SUBROUTINE momentsc(n, p, fval)
  USE commonvars

  IMPLICIT NONE
  INCLUDE 'mpif.h'

  INTEGER, INTENT(IN)     :: n
  REAL(8), INTENT(IN)     :: p(n)
  REAL(8), INTENT(OUT)    :: fval
  INTEGER                 :: i, tt, imun, ii, npoly
  REAL(8)                 :: temp(Nmom,Nmomdata), diffmom(Nmom,1), qcons_w_riv, t1, t2, sumbetadf, uu, pfine_temp, pfine_temp1, cdf_fine(Nfines), true_cdf_fine(Nfines), frac_fine_temp, &
                             prob_a_temp, prob_a_temp1, cdf_a(Nab), true_cdf_a(Nab), a_discr_temp, beta, betap, sigma_ap

  EXTERNAL                :: utilfun, emaxfun, simulatemun

  CALL cpu_time(t1)

  iterations=iterations+1
  IF(myrank == 0) write(*,*) 'iteration', iterations

  ii = 0
  
  ii = ii + 1
  !theta is the relative value of public consumption; we switched to a model with only one type in terms of relative value
  IF (ii <= nparam) THEN
     theta = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  ELSE
     theta = 0.1d0
  END IF
  
  ii = ii + 1
  ! Electoral appeal parameter
  IF (ii <= nparam) THEN
     alpha_ea = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  ELSE
     alpha_ea = 0.5d0
  END IF

  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     alpha_st = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  ELSE
     alpha_st = -0.5d0
  END IF
  
  ii = ii + 1  
  !probability of a high electoral appeal mayor
  IF (ii <= nparam) THEN
     pappeal(2) = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  ELSE
     pappeal(2) = 0.5d0
  END IF

  ii = ii + 1
  IF (ii <= nparam) THEN
     cost_run_frac = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  ELSE
     cost_run_frac = 0.1d0
  END IF

  ii = ii + 1
  IF (ii <= nparam) THEN
     upower = cost_run_frac + factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  ELSE
     upower = cost_run_frac + 1d0
  END IF

  ii = ii + 1
  IF (ii <= nparam) THEN
     cost_run_frac_st = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  ELSE
     cost_run_frac_st = 0.05d0
  END IF

  ii = ii + 1
  IF (ii <= nparam) THEN
     sigma_steal = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  ELSE
     sigma_steal = 0.01d0
  END IF
  
  !params - no bounds
  ! ELECTORAL RULE WITH CORRECT AGE VARIABLE: 9 params, but different vars:
  !Constant in pmayor
  alpha(1) = -5.1205603d0 !-5.2217137d0 !(ubp(1)*EXP(p(1))+lbp(1))/(1d0+EXP(p(1)))
  !Parameter on the audit dummy in pmayor
  alpha(2) = .3281086d0 !(ubp(2)*EXP(p(2))+lbp(2))/(1d0+EXP(p(2)))
  !Parameter on stealing multiplied by the audit dummy in pmayor
  alpha(3) = -.5535916d0 !(ubp(3)*EXP(p(3))+lbp(3))/(1d0+EXP(p(3)))
  !Parameter on log public consumption (public gdp) in pmayor
  alpha(4) = .5534424d0 !(ubp(4)*EXP(p(4))+lbp(4))/(1d0+EXP(p(4)))
  !Parameter on age in pmayor
  alpha(5) = -.1276038d0 !(ubp(5)*EXP(p(5))+lbp(5))/(1d0+EXP(p(5)))
  !Parameter on high electoral appeal in pmayor
  alpha(6) = 4.8565301d0 !(ubp(6)*EXP(p(6))+lbp(6))/(1d0+EXP(p(6)))
  !Parameter on dummy for large municipality
  alpha(7) = .6110852d0 !(ubp(7)*EXP(p(7))+lbp(7))/(1d0+EXP(p(7)))
  !The parameter on log population size is equal to alpha(6)*eta
  !Standard Deviation in pmayor
  !THE STANDARD DEVIATION MUST BE NORMALIZED TO 1
  alpha(8) = 1d0
  
  !probability of a low electoral appeal mayor
  pappeal(1) = 1d0 - pappeal(2)
  !low electoral appeal
  ptype(1) = pappeal(1)
  !high electoral appeal
  ptype(2) = pappeal(2)

!!$. ml model lf electoral_rule_lf (xb: reeleito00_04_b = audit d_fr_stolen_audit lngdp_01_04 lngdp_priv_01_04 mayor_age relcon2 ln_pop2001 PesResAlfaRate pop_d3) (ln_const:) (tr_p:), technique(bhhh)
!!$-----------------------------------------------------------------------------------
!!$                  |                 OPG
!!$  reeleito00_04_b |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
!!$------------------+----------------------------------------------------------------
!!$xb                |
!!$            audit |   .3281086   .2386803     1.37   0.169    -.1396961    .7959133
!!$d_fr_stolen_audit |  -.5535916   .2743165    -2.02   0.044    -1.091242   -.0159411
!!$      lngdp_01_04 |   .5534424   .1685504     3.28   0.001     .2230897    .8837951
!!$ lngdp_priv_01_04 |  -.3358938   .1126018    -2.98   0.003    -.5565893   -.1151983
!!$        mayor_age |  -.1276038   .0184315    -6.92   0.000    -.1637289   -.0914787
!!$          relcon2 |     .32518   .0426875     7.62   0.000      .241514    .4088459
!!$       ln_pop2001 |  -.4689685   .0998538    -4.70   0.000    -.6646783   -.2732586
!!$   PesResAlfaRate |  -.0001188    .000046    -2.58   0.010     -.000209   -.0000287
!!$           pop_d3 |   .6110852   .1623053     3.77   0.000     .2929726    .9291978
!!$            _cons |  -1.762454   1.733389    -1.02   0.309    -5.159835    1.634927
!!$------------------+----------------------------------------------------------------
!!$ln_const          |
!!$            _cons |   1.580324   .3233808     4.89   0.000     .9465094    2.214139
!!$------------------+----------------------------------------------------------------
!!$tr_p              |
!!$            _cons |   2.186069   .3707693     5.90   0.000     1.459375    2.912764
!!$-----------------------------------------------------------------------------------
!!$. di gamma_0_1
!!$-1.762454
!!$
!!$. di gamma_0_2
!!$3.0940761
!!$
!!$. di gamma_0_1_adj
!!$-5.1205603
!!$
!!$. di gamma_0_2_adj
!!$-.26403018
!!$
!!$. di value_ht
!!$4.8565301
!!$
!!$. di p_gamma
!!$.89899152
!!$
!!$. di one_p_gamma
!!$.10100848
!!$
!!$. di gamma_1
!!$.55344242
!!$
!!$. di eta
!!$.84736632

  !We use a log-normal distribution for the fine
  !mean of the fine distribution
  mu_fine = 0.0942693d0 ! (ubp(1)*EXP(p(1))+lbp(1))/(1d0+EXP(p(1)))                                                                      
  !standard deviation of the fine distribution
  sigma_fine = 0.2837051d0 ! (ubp(2)*EXP(p(2))+lbp(2))/(1d0+EXP(p(2)))

  IF (flag_mean_fine == 1) sigma_fine = 0d0 

  !The following calculations produce the discrete approximation for the probability of drawing one of the discrite values, frac_fine(:), and one of the corresponding discrete values, pfine(:)
  !Given the parameters of the fine distribution we derive the support points that gives the best approximation (part B in Kennan)
  DO i = 1, Nfines
     ! cdf values at the fixed support points
     true_cdf_fine(i) = (2d0*DBLE(i)-1)/(2d0*DBLE(Nfines))
     !We compute the inverse of the uniform at true_cdf_fine(i); we consider a uniform defined on [0,mu_fine+sigma_fine*E[stealing in data]]
     CALL normal_01_cdf_inv(true_cdf_fine(i),frac_fine_temp)
     frac_fine(i) = EXP(sigma_fine*frac_fine_temp + mu_fine)
  END DO

  !We compute the probabilities at the chosen points that best approximate the log-normal distribution (part A in Kennan)
  CALL normal_01_cdf((LOG(frac_fine(1)) - mu_fine)/sigma_fine,pfine_temp)
  DO i = 2, Nfines
     CALL normal_01_cdf((LOG(frac_fine(i)) - mu_fine)/sigma_fine,pfine_temp1)
     cdf_fine(i-1) = (pfine_temp + pfine_temp1)/2d0
     pfine_temp = pfine_temp1
  END DO
  cdf_fine(Nfines) = 1d0

  !We compute the probability of the discrete fine values in our grid
  pfine(1) = cdf_fine(1)
  DO i = 2, Nfines
     pfine(i) = cdf_fine(i) - cdf_fine(i-1)
  END DO

  IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf_fine', true_cdf_fine(:)
  IF(myrank == 0) write(*,'(a)') ''
  IF(myrank == 0) write(*,'(a,3f20.10)') 'frac_fine', frac_fine(:)
  IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf_fine', cdf_fine(:)
  IF(myrank == 0) write(*,'(a,3f20.10)') 'pfine', pfine(:)
  
  !The ability of mayors is drawn from LOG-N(mu_a,sigma_a)
  !Mean of the ability distribution, we normalize the constant in the production function to 0 and assign the entire constant in the IV estimation to mu_a
  mu_a = -0.2941902d0 !-0.32856629d0 !(ubp(1)*EXP(p(1))+lbp(1))/(1d0+EXP(p(1))) !-.53632455d0
  !standard deviation of the ability distribution
  sigma_a = 0.39916454d0 !0.44185212d0 ! .35204222d0 !.3438509d0 !.4730158d0 ! .48566405d0
  !standard deviation of the ability distribution in the second term
  sigma_ap = 0.46006851d0 !0.48344978d0
  !Mean of the ability distribution in the second term = mu_a + beta' - beta, where beta' and beta are the estimated constants in the production functions for the first and second term
  beta = -.2941902d0 !-0.32856629d0
  betap = -.2829773d0 !-0.30641002d0

  !parameters of the production function for public consumption                                         
  !we use the same constant in both terms
  !power parameter for public inputs
  alphaq(1) = 0.1526144d0 !0.1518577d0
  !power parameter for private inputs
  alphaq(2) = 0.429166d0 !0.3725852d0
  !linear term in education
  alphaq(3) = 0.0430096d0 !0.0340758d0
  !We normalize the constant to be equal to 0
  alphaq(4) = beta - mu_a
  !flexible polynomial in population
  alphaq(5) = 6.66d-06 
  alphaq(6) = -2.26d-10
  alphaq(7) = 1.46d-15
  alphaq(8) = -3.75d-21
  alphaq(9) = 3.16d-27

  !Polynomial of order n in population
  Npoly = Nparq - 4
  polypop = 0d0
  DO imun = 1, Nmun
     DO i = 1, Npoly
        polypop(imun) = alphaq(4+i)*popsize(imun)**DBLE(i)
        !if (myrank ==0) write(*,'(a,2i4,4f20.10)') 'polypop', imun, i, polypop(imun), alphaq(4+i), popsize(imun), popsize(imun)**DBLE(i)
     END DO
  END DO

  !Ea and Eap are the expectation of ability in the first and second term (ability is log-normally distributed)
!!$  Ea = EXP(mu_a+sigma_a**2d0/2d0)
!!$  Eap = EXP(mu_ap+sigma_ap**2d0/2d0)

!!$. ivregress gmm lngdp_pc_01_04 lindustry_index pref_edu PesResAlfaRate* pib2001* d_2nd_term (lpinputs_pc_01_04 = fpm_d1_200* fpm_d2_200* fpm_d3_200*  fpm_d4_200*  fpm_d5_20
!!$> 0* fpm_d6_200* fpm_d7_200* fpm_d8_200* fpm_d9_200* fpm_d10_200* fpm_d11_200* fpm_d12_200* fpm_d1_pop_200* employment_rate2000 wage_priv nfirms ln_pop_01_04) if nterms >= 
!!$> 1, rob
!!$ Instrumental variables (GMM) regression           Number of obs   =        474
!!$                                                  Wald chi2(8)    =    7598.90
!!$                                                  Prob > chi2     =     0.0000
!!$                                                  R-squared       =     0.6192
!!$GMM weight matrix: Robust                         Root MSE        =     .46004
!!$
!!$------------------------------------------------------------------------------------
!!$                   |               Robust
!!$    lngdp_pc_01_04 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
!!$-------------------+----------------------------------------------------------------
!!$ lpinputs_pc_01_04 |   .1518577   .0379653     4.00   0.000     .0774472    .2262683
!!$   lindustry_index |   .3725852    .082993     4.49   0.000     .2099218    .5352485
!!$          pref_edu |   .0340758   .0377231     0.90   0.366    -.0398601    .1080118
!!$    PesResAlfaRate |   .0003879   .0000213    18.23   0.000     .0003462    .0004296
!!$PesResAlfaRate_2nd |  -.0000556   .0000235    -2.36   0.018    -.0001017   -9.51e-06
!!$           pib2001 |   1.32e-07   9.36e-08     1.41   0.159    -5.15e-08    3.15e-07
!!$       pib2001_2nd |   6.25e-07   1.87e-07     3.35   0.001     2.59e-07    9.91e-07
!!$        d_2nd_term |   .3983394   .1729372     2.30   0.021     .0593888    .7372901
!!$             _cons |  -3.391517   .2256809   -15.03   0.000    -3.833843    -2.94919
!!$------------------------------------------------------------------------------------
!!$
!!$. di beta_1
!!$-.32856629
!!$. di beta_2
!!$-.30641002
!!$. di std_dev_all_t1
!!$.44185212
!!$. di std_dev_all_t2
!!$.48344978

  !The following calculations produce the discrete approximation for the probability of drawing one of the discrite values, prob_a(:), and one of the corresponding discrete values, a_discr(:)
  !Given the parameters of the ability distribution we derive the support points that gives the best approximation (part B in Kennan)
  DO i = 1, Nab
     ! cdf values at the fixed support points
     true_cdf_a(i) = (2d0*DBLE(i)-1)/(2d0*DBLE(Nab))
     !We compute the inverse of the standard normal at true_cdf_a(i)
     CALL normal_01_cdf_inv(true_cdf_a(i), a_discr_temp)
     a_discr(i) = sigma_a*a_discr_temp + mu_a
  END DO
  
  !We compute the probabilities at the chosen points that best approximate the log-normal distribution (part A in Kennan)
  CALL normal_01_cdf((a_discr(1) - mu_a)/sigma_a,prob_a_temp)
  DO i = 2, Nab
     CALL normal_01_cdf((a_discr(i) - mu_a)/sigma_a,prob_a_temp1)
     cdf_a(i-1) = (prob_a_temp + prob_a_temp1)/2d0
     prob_a_temp = prob_a_temp1
  END DO
  cdf_a(Nab) = 1d0

  !We compute the probability of the discrete fine values in our grid
  pability(1) = cdf_a(1)
  DO i = 2, Nab
     pability(i) = cdf_a(i) - cdf_a(i-1)
  END DO

  !We take the exp of a_discr(i) to make it a log-normal
  a_discr(:) = EXP(a_discr(:))

  IF(myrank == 0) write(*,'(a)') ''
  IF(myrank == 0) write(*,'(a,3f20.10)') 'true cdf ability', true_cdf_a(:)
  IF(myrank == 0) write(*,'(a,3f20.10)') 'discrete ability', a_discr(:)
  IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf ability', cdf_a(:)
  IF(myrank == 0) write(*,'(a,3f20.10)') 'prob discrete ability', pability(:)
    
  !parameters of the wage function for past mayors
  !constant
  alphaw(1) = -3.086317d0
  !education
  alphaw(2) = 4.096486d0
  !age
  alphaw(3) = 1.035222d0
  !age^2
  alphaw(4) = -.0975276d0
  !audit*ln(popsize)
  alphaw(5) = 0d0
  !audit*d_steal*ln(popsize)
  alphaw(6) = 0d0
  !medium size municipality
  alphaw(7) = 1.144047d0
  !large size municipality
  alphaw(8) = 3.215137d0
  !standard deviation
  sigmaw = 6.335244d0 !6.464272d0

!!$  ! If stealing has no effect on future wages, we don't need to run the loops over ipast in emaxfun and emaxfuntemp
!!$  IF ((alphaw(5)>-0.0001d0 .AND. alphaw(5)<0.0001d0) .AND. (alphaw(6)>-0.0001d0 .AND. alphaw(6)<0.0001d0)) Npast = 1
!!$  ! If we simulate the policy with a tax on future wages, we need Npast = 3
!!$  IF (flag_tax == 1) Npast = 1

!!$ reg wage_mean_past_norm  pref_edu age age2 pop_d2 pop_d3 if wage_mean_past_norm >.5668408  & wage_mean_past_norm <  10000 & age<99, rob
!!$
!!$Linear regression                               Number of obs     =      3,725
!!$                                                F(5, 3719)        =     143.06
!!$                                                Prob > F          =     0.0000
!!$                                                R-squared         =     0.1392
!!$                                                Root MSE          =     6.3404
!!$
!!$------------------------------------------------------------------------------
!!$             |               Robust
!!$wage_mean_~m |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
!!$-------------+----------------------------------------------------------------
!!$    pref_edu |   4.096486   .1794472    22.83   0.000     3.744662    4.448311
!!$         age |   1.035222   .1830252     5.66   0.000     .6763825    1.394062
!!$        age2 |  -.0975276   .0182942    -5.33   0.000    -.1333952     -.06166
!!$      pop_d2 |   1.144047    .211215     5.42   0.000     .7299387    1.558156
!!$      pop_d3 |   3.215137   .3393853     9.47   0.000     2.549737    3.880536
!!$       _cons |  -3.086317   .5092068    -6.06   0.000    -4.084668   -2.087965
!!$------------------------------------------------------------------------------

!!$  !constant
!!$  alphaw(1) = -2.165371d0
!!$  !education
!!$  alphaw(2) = .6182883d0
!!$  !age
!!$  alphaw(3) = .3021505d0
!!$  !age^2
!!$  alphaw(4) = -.0088004d0
!!$  !audit*ln(popsize)
!!$  alphaw(5) = 0d0
!!$  !audit*d_steal*ln(popsize)
!!$  alphaw(6) = 0d0
!!$  !popsize
!!$  alphaw(7) = .1967507d0
!!$  !popsize**2
!!$  alphaw(8) = .5046361d0
!!$  !standard deviation
!!$  sigmaw = .8073561d0 !6.464272d0

!!$. reg lwage_mean_past_norm  pref_edu age age2 missing_age  pop_d2-pop_d3, rob
!!$
!!$Linear regression                               Number of obs     =      5,334
!!$                                                F(6, 5327)        =     250.45
!!$                                                Prob > F          =     0.0000
!!$                                                R-squared         =     0.2141
!!$                                                Root MSE          =     .80789
!!$
!!$------------------------------------------------------------------------------
!!$             |               Robust
!!$lwage_mean.. |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
!!$-------------+----------------------------------------------------------------
!!$    pref_edu |   .6182883   .0237276    26.06   0.000     .5717724    .6648041
!!$         age |   .3021505    .066254     4.56   0.000     .1722655    .4320354
!!$        age2 |  -.0088004   .0021045    -4.18   0.000    -.0129261   -.0046747
!!$ missing_age |   58.62389   14.62151     4.01   0.000     29.95974    87.28804
!!$      pop_d2 |   .1967507   .0250966     7.84   0.000     .1475511    .2459502
!!$      pop_d3 |   .5046361   .0319755    15.78   0.000     .4419511    .5673211
!!$       _cons |  -2.165371   .5174155    -4.18   0.000    -3.179717   -1.151025
!!$------------------------------------------------------------------------------

  !Standard deviation of the private input shocks
  sigma_zzpr = 0d0 !(ubp(20)*EXP(p(20))+lbp(20))/(1d0+EXP(p(20)))

  !Standard deviation of the shocks to the decision to run
  sigma_run = 0.1d0 ! 1d0 !0.1d0 ! (ubp(19)*EXP(p(19))+lbp(19))/(1d0+EXP(p(19)))
      
  !Degree of Rivalry Parameter
!!$  eta = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
  eta = 1d0

  ! WE WILL ASSUME LOG PREFERENCES
  !preference parameters                                                                                
  delta = 1d0 ! (ubp(6)*EXP(p(6))+lbp(6))/(1d0+EXP(p(6))) if we use log utility for public consumption, otherwise we have to estimate the parameter 
  !gamma = 1 if we use log utility for private consumption, otherwise we have to estimate the parameter 
  gamma = 2d0

  ! WE SHOULD GET RID OF THE MEASURMENT ERRORS
  !Mean of the measurement error in stealing; a positive mean means that it is more likely that a positive clearical mistake (leading to corruption) is made than negative (leading to negative corruption)
  !As a percentage of funds (see simulatemun)
  mu_steal = 0d0 !(ubp(22)*EXP(p(22))+lbp(22))/(1d0+EXP(p(22)))
  
  ! I DON'T UNDERSTAND WHAT THIS PARAMETER DOES
  !We transform the level of public consumption to take into account the degree of rivalry
  !This is also the variable we output for public consumption
  !If we divide by (popsizeg)**eta, the utility at subsistence level goes to -infinity 
  qcons_w_riv = subqcons + small_no

  subutil = 0d0
  CALL utilfun(subpcons,qcons_w_riv,theta,uu)
  IF (myrank==0) WRITE(*,*) 'uu,subpcons,subqcons,theta', uu,subpcons,qcons_w_riv,theta
  sumbetadf = 0d0
  DO tt = 1, Tend
     sumbetadf = sumbetadf + betadf**(tt-1)
     subutil(Tend - tt + 1) = uu*sumbetadf
  END DO

!!$  subutil = (/ -9146.35964062422, -8312.61187986127, -7461.84885867459, -6593.72332685144, -5707.88094744007, -4803.96015212234, -3881.59199363487, -2940.39999517826, -1979.99999675314, -999.999998360174 /)

  IF (myrank==0) WRITE(*,*) 'subutil', subutil(:)

  IF (myrank == 0) THEN
     write(*,*) ''
     write(*,*) 'param', p(:)
     write(*,'(a,19f12.8)') 'pp', (/ alpha_ea, alpha_st, mu_fine, sigma_fine, ptype(1), alphaq(1:9), theta, sigma_zzpr, cost_run_frac, upower, eta /) 
     write(*,*) 'electoral appeal', alpha_ea, alpha_st
     write(*,*) 'mu_fine, sigma_fine', mu_fine, sigma_fine
     write(*,*) 'prob low electoral appeal', ptype(1)
     write(*,*) 'prod fn: pu*ab, pr*ab, educ, const', alphaq(1:4)
     write(*,*) 'prod fn: poly pop', alphaq(5:Nparq)
     write(*,*) 'theta', theta
     write(*,*) 'wages past mayors', alphaw, sigmaw
     write(*,*) 'delta, gamma', delta, gamma
     write(*,*) 'sigma_run, sigma_zzpr', sigma_run, sigma_zzpr
     write(*,*) 'mu_steal, sigma_steal', mu_steal, sigma_steal
     write(*,*) 'cost_run_frac, cost_run_frac_st, upower, upower-cost_run_frac', cost_run_frac, cost_run_frac_st, upower, upower-cost_run_frac
     write(*,*) 'high appeal param, prob of low and high appeal', alpha(7), pappeal(:)
     write(*,*) 'mu_a, sigma_a', mu_a, sigma_a
     write(*,*) 'beta, betap', beta, betap
     write(*,*) 'eta', eta
     write(*,'(a,8f18.14)') 'param inputs', alpha_ea, alpha_st, pappeal(2), theta, cost_run_frac, upower-cost_run_frac, cost_run_frac_st, sigma_steal
  END IF
  
  CALL emaxfun()
  
  if (myrank == 0) write(*,*) 'here 1'
  
  CALL simulatemun()

  DEALLOCATE(Pr_i_win)

  if (myrank == 0) write(*,*) 'here 3'
  
  !----------- COMPUTE CRITERION FUNCTION ------------!

  ! Compute the difference in moments
  diffmom(1:Nmom,1) = simmom(1:Nmom,1) - datamom(1:Nmom,1)

  DO i=1,Nmom
     temp(i,1)=SUM(diffmom(:,1)*INV_V(:,i))
     if (myrank == 0) write(*,'(A,1I3,x,A,1F20.8,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6)') 'i',i,'dfval',diffmom(i,1)*INV_V(i,i)*diffmom(i,1),'simmom',simmom(i,1),'datamom',datamom(i,1),'diff',diffmom(i,1),'V',V(i,i),'weight',INV_V(i,i)
  END DO
  if (myrank == 0) write(*,*) ''

  fval=SUM(temp(:,1)*diffmom(:,1))

  CALL cpu_time(t2)
  IF (myrank == 0) THEN
     write(*,*) 'simmom', simmom
     write(*,*) 'datamom', datamom
     write(*,*) 'fval', fval
     write(*,*) 'time for one set of parameters', t2-t1
  END IF

  IF (flag_simulate_data == 1) THEN

     STOP
     CALL MPI_FINALIZE(ierr)

  END IF

END SUBROUTINE momentsc
